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ABSTRACT 

We compare the measured angular cross-correlation between the Fermi-LAT 7-ray sky and catalogues of 
extra-galactic objects with the expected signal induced by weakly interacting massive particle (WIMP) dark 
matter (DM). We include a detailed description of the contribution of astrophysical 7-ray emitters such as 
blazars, misaligned AGN and star forming galaxies, and perform a global ht to the measured cross-correlation. 

Live catalogues are considered; SDSS-DR6 quasars, 2MASS galaxies, NVSS radio galaxies, SDSS-DR8 Lu¬ 
minous Red Galaxies and SDSS-DR8 main galaxy sample. To model the cross-correlation signal we use the 
halo occupation distribution formalism to estimate the number of galaxies of a given catalogue in DM halos 
and their spatial correlation properties. We discuss uncertainties in the predicted cross-correlation signal aris¬ 
ing from the DM clustering and WIMP microscopic properties, which set the DM 7-ray emission. The use 
of different catalogues probing objects at different redshifts reduces significantly, though not completely, the 
degeneracy among the different 7-ray components. We hnd that the presence of a significant WIMP DM signal 
is allowed by the data but not significantly preferred by the fit, although this is mainly due to a degeneracy 
with the misaligned AGN component. With modest substructure boost, the sensitivity of this method excludes 
thermal annihilation cross sections at 95% level, for WIMP masses up to few tens of GeV. Constraining the 
low-redshift properties of astrophysical populations with future data will further improve the sensitivity to DM. 

Subject headings: cosmology: theory - cosmology: observations - cosmology; dark matter - cosmology: 
large-scale structure of the universe - gamma rays: diffuse backgrounds 


1. INTRODUCTION 

The last few years have seen a tremendous improvement 
in our understanding of the 7-ray sky, mostly thanks to the 
observations performed by the Large Area Telescope (EAT) 
on board of the Fermi satellite (Atwood et al. 2009). Among 
the main issues that have been investigated, an important one 
is the understanding of the origin of the Isotropic Gamma- 
Ray Background (IGRB) (Ackermann et al. 2015a; Eornasa 
(& Sanchez-Conde 2015), i.e. the fraction of the extra-galactic 
7-ray background (EGB) that has not been resolved into in¬ 
dividual sources. The nature of the extragalactic emission 
is a recurrent issue which arises each time a new observa¬ 
tional window of the electromagnetic spectrum is opened on 
the Universe. A good example is the quest for the origin of 
the soft X-ray background, with the important difference that 
the latter has now been largely resolved (see, e.g., Hickox & 
Markevitch (2007)) whereas a significant fraction of the 7- 
ray flux is still diffuse, leaving large room for potential new 
discoveries. 

The interest in the IGRB also stems from the consideration 
that the 7-ray band is a potential “golden channel" for the in¬ 
direct detection of particle Dark Matter (DM). In fact, among 
the conventional astrophysical sources that contribute to the 
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IGRB there is the possibility that a characteristic signal from 
DM annihilation or decay may also be present. After its first 
detection and early attempts to shed light on the origin of the 
IGRB (see e.g. Kraushaar et al. (1972); Eichtel et al. (1973); 
Mayer-Hasselwander et al. (1982); Padovani et al. 0993); 
Stecker & Salamon (1996); Sreekumar et al. (1998); Keshet 
et al. (2004); Strong et al. (2004)), a significant step forward 
has recently been possible thanks to Fermi-LAT that is resolv¬ 
ing an ever growing number of sources (Nolan et al. 2012; 
Acero et al. 2015; Ackermann et al. 2015b), most of which 
have been identified as blazars, almost equally split into Flat 
Spectrum Radio Quasars (FSRQs) and BE Lacs sub-classes. 

The properties of the resolved sources can be used to ex¬ 
trapolate their contribution to the IGRB (Ajello et al. 2012, 
2014; Di Mauro et al. 2014c,b; Ajello et al. 2015). These 
population studies suggest that unresolved blazars account 
for only about 20% of the unresolved IGRB integrated above 
100 MeV, while they can be the dominant component above 
few GeV. The remaining IGRB fraction is thought to be con¬ 
tributed by star-forming galaxies (SFGs) and misaligned AGN 
(mAGN), two types of sources that can contribute 10-50% 
each to the extragalactic 7-ray emission (Inoue 2011; Acker¬ 
mann et al. 2012; Di Mauro et al. 2014a). The contribution 
from additional potential sources like the millisecond pulsars 
located in our Galaxy at high Galactic latitudes turned out to 
be small (Siegal-Gaskins et al. 2011; Galore et al. 2014). The 
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contribution from known astrophysical sources to the IGRB 
has thus significant uncertainties and leaves room for an addi¬ 
tional contribution by more exotic sources like DM. 

Additional constraints on the origin of the IGRB can be ob¬ 
tained by analyzing the angular correlation properties of its 
fluctuations. Refs. Ackermann et al. (2012); Cuoco et al. 
(2012); Ackermann et al. (2012) have confirmed the conclu¬ 
sions derived from the IGRB mean-intensity and source popu¬ 
lations studies: a blazar population that contributes, at energy 
below ^10 GeV, about 20% of the unresolved IGRB can ac¬ 
count for the whole measured angular power, thus providing 
an independent confirmation that the IGRB is not dominated 
by emission from blazars in the low energy part of the spec¬ 
trum. Constraints on the DM contribution have been derived 
in Ando & Komatsu (2013); Gomez-Vargas et al. (2014); For- 
nasa et al. (2013). 

The accuracy in the analysis of the IGRB and its fluctua¬ 
tions is limited by the presence of Galactic foregrounds and 
bright sources. If incorrectly subtracted they can induce spu¬ 
rious contributions to both the mean IGRB intensity and its 
anisotropies. An effective way of dealing with this problem 
and Alter out contaminations is to cross-correlate the IGRB 
with maps of sources (observed in other wavelengths or by 
other means) that trace the same structures where the actual 
IGRB sources reside but do not correlate with Galactic fore¬ 
grounds. Basically, all catalogs of extragalactic objects at any 
redshift satisfy these conditions. In the framework of the 
IGRB investigation, the cross-correlation strategy has been 
proposed in Cuoco et al. (2008); Ando & Pavlidou (2009) and 
recently revisited in Ando (2014); Ando et al. (2014). The 
measurement was pioneered in Xia et al. (2011) using the first 
21 months of Fermi data. In that case, no statistically signifi¬ 
cant signal was observed. The analysis has then been recently 
updated using 60-months Fermi maps (Xia et al. 2015). This 
time a significant (more than 3.5 cr C.L.) cross-correlation sig¬ 
nal has been detected. The signal is present on angular scales 
smaller than 1 ° in the cross-correlation between the diffuse 
7 -ray emission cleaned by the Galactic foregrounds and four 
types of Large Scale Structure (LSS) tracers: radio galax¬ 
ies in the NRAO VLA Sky Survey (NVSS) (Blake & Wall 
2002), near infra-red selected galaxies in the Two Micron All 
Sky Survey (2MASS) (Jarrett et al. 2000) , optically selected 
galaxies in the Sloan Digital Sky Survey (SDSS)-DR 8 cat¬ 
alog (Aihara et al. 2011) and quasi stellar objects (QSO) in 
the SDSS-DR 6 catalogs (QSO-DR 6 ) (Richards et al. 2009). 
No significant correlation was observed with Luminous Red 
Galaxies (LRG) also from SDSS. The analysis further con¬ 
firms that blazars provide a minor contribution < 20 % to the 
IGRB as found in the IGRB mean-intensity and source pop¬ 
ulations studies, while a mixture of SFGs and mAGNs can in 
principle contribute to the majority of the IGRB. 

A promising, possibly more effective in the context of DM, 
way to apply the cross-correlation technique is to use weak 
gravitational lensing maps (cosmic shear) instead of cata¬ 
logs of LSS tracers (Camera et al. 2013; Fomengo & Regis 
2014; Fomengo et al. 2014; Camera et al. 2014; Shirasaki 
et al. 2014). This alternative approach, originally proposed in 
(Camera et al. 2013), has the advantage of probing directly the 
matter distribution, therefore avoiding the so-called ‘biasing’ 
issue, i.e., the fact that the mapping between the spatial distri¬ 
bution of extragalactic sources and that of the underlying mass 
density held is ill-known, and needs to be modelled. Cross¬ 
correlation of 7 -rays with cosmic shear will become available 
in the next years with the release of cosmic-shear maps from 


wide area surveys, like, e.g., the Dark Energy Survey (DES) 
(The Dark Energy Survey Collaboration 2005) and, in the next 
decade, by the satellite-based Euclid survey (Laureijs 2009). 
Einally, a similar technique, based on the cross-correlation of 
7 -rays with Cosmic Microwave Background (CMB) lensing 
maps, has been recently adopted in Eornengo et al. (2014), 
where an evidence of 3.2a has been reported providing a fur¬ 
ther direct evidence of the extragalactic origin of the IGRB, 
and of a subdominant role of Galactic sources. 

In this paper we investigate the implications of the recent 
measurement of a cross-correlation between 7 -rays and LSS 
tracers by Xia et al. (2015) for both the DM and the main as¬ 
trophysical contributors to the IGRB. This work builds upon 
the results obtained by Regis et al. (2015) in which we con¬ 
centrated on the low-redshift 2MASS catalog as a tracer of 
the LSS in the local Universe, and we have assumed that 
DM-induced 7 -rays provide the dominant source of cross¬ 
correlation for such a low redshift observations. That ap¬ 
proach has been motivated by the fact that the DM contri¬ 
bution to the cross-correlation is dominated by 7 -rays emis¬ 
sion at low-redshift (see e.g. Eornengo & Regis (2014) or 
Appendix A), which is where 2MASS galaxies mostly reside. 
In that analysis we found that the observed cross-correlation 
signal can indeed be explained by a DM emission, while its 
contribution to the total mean intensity is significantly below 
the IGRB intensity measured by Fermi. This implies that the 
cross-correlation technique can be a powerful probe of the 
particle nature of DM, even when the DM contribution to the 
IGRB is subdominant, which is what we expect in a realis¬ 
tic scenario. In (Regis et al. 2015) we found that the cross¬ 
correlation signal can be explained by a DM particle with 
mass in the tens to hundreds GeV range (depending on the 
7 -rays production channel) and, once the uncertainties in the 
DM distribution modeling is properly accounted for, a “ther¬ 
mal” value for the annihilation cross section ((ctov) = 3 x 10“^® 
cm^ s“^), which is the most appealing case for a weakly inter¬ 
acting massive particle (WIMP) DM. Erom the same analysis 
we have obtained upper bounds on the DM annihilation cross 
section and decay rate that turn out to be quite competitive 
with those obtained with different techniques, based either on 
local (Galactic halo, dwarf galaxies) or extra-galactic 7 -ray 
emission. We point out that those constraints are conservative 
precisely because the DM is assumed to be the only source of 
the 7 -ray signal. 

In this follow-up paper we extend the study of Regis et al. 
(2015) to the inclusion of astrophysical 7 -ray emitters, and 
to the whole set of LSS-tracers catalogs. As it will be dis¬ 
cussed in the next Sections, the redshift distributions of the 
7 -ray signal is a fingerprint that characterises the contribution 
of different astrophysical sources and of the DM. Eor this rea¬ 
son, the possibility to use catalogues of objects whose distri¬ 
butions peak at different redshifts is an effective way to extract 
the information encoded in the 7 -rays maps and remove de¬ 
generacies. To this aim, in addition to DM, we account here 
for contributions from blazars (of both BL Lac and ESRQs 
types), SEGs and mAGNs. Consequently, we do not limit 
our cross correlation study to the 2MASS catalog but con¬ 
sider NVSS, SDSS-DR 6 QSO, LRGs and SDSS-DR 8 Main 
Galaxy samples as well. The approach will be similar to Xia 
et al. (2015) where, however, only astrophysical sources have 
been fitted to the observed correlations. Here, beside includ¬ 
ing DM in the fit, we will use an improved description of 
the cross-correlation modeling between astrophysical sources 
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and LSS tracers based on the Halo Occupation Distribution 
(HOD) formalism. As for the DM, we shall use the halo 
model to trace its spatial distribution and predict its cross¬ 
correlation with LSS (see e.g., Cooray & Sheth (2002); Ando 
& Komatsu (2013); Fornengo & Regis (2014)). 

The paper is organized as follows. In Section 2 we present 
the theoretical estimate of the angular cross-correlation func¬ 
tion and angular power spectra. Section 3 describes the statis¬ 
tical techniques employed in the determination of the param¬ 
eters of the 7 -rays emitters (DM and astrophysical sources) 
from the measured cross-correlation reported in Xia et al. 
(2015). Section 4 then shows our results, and finally Sec¬ 
tion 5 summarizes our conclusions. The technical aspects of 
our theoretical modeling are presented in a set of three Ap¬ 
pendices. Appendix A introduces the modeling of the win¬ 
dow functions of DM and astrophysical 7 -rays sources and 
of catalogs of LSS tracers. Appendix B discusses the HOD 
of galaxies for the various catalogs. Appendix C describes 
the derivation of three-dimensional (3D) power spectra (PS). 
These are the ingredients used in Section 2. 

In this work we assume a fiducial flat ACDM model with 
the cosmological parameters derived by the Planck Collab¬ 
oration in Planck Collaboration et al. (2015): matter density 
parameter = 0.31, baryon density parameter = 0 . 022 , 
reduced Hubble constant h = 0.68, rms matter fluctuations in 
a comoving sphere of 8 Mpc erg = 0.83 and spectral index of 
primordial scalar perturbations = 0.96. 

2 . FORMALISM 

To quantify the cross-correlation between 7 -ray sources and 
the LSS tracers in the various catalogues, we consider both 
the 2-point angular cross-correlation function (CCF) and its 
Legendre transform, i.e. the cross angular power spectrum 
(CAPS). In the Limber approximation (Limber 1953; Kaiser 
1992, 1998), the CAPS can be obtained by integrating the 3D- 
PS of cross-correlation P^g(k,z): 

J X 

where x(z) denotes the radial comoving distance, W(x) is the 
so-called window function that characterizes the distribution 
of objects and 7 -ray emitters along the line of sight, k is the 
modulus of the wavenumber and I is the multipole. We relate 
the cosmological redshift z to the radial comoving distances 
X through the differential relation, valid in a flat cosmology, 
dx = cdz/H{z), where H{z) is the expansion rate of the Uni¬ 
verse. 

The indices 7 and g denote 7 -ray emitters and extragalactic 
objects in different catalogs, respectively. We consider five 
types of 7 -ray sources: three different flavours of AGNs (BL 
Lacs, FSQRs, mAGN), SFGs and DM. We will consider both 
the case of annihilating and decaying DM particles. For the 
LSS tracers, we consider five different catalogues: quasars in 
SDSS-DR 6 , 2MASS galaxies, NVSS radio sources, SDSS- 
DR 8 Luminous Red Galaxies and SDSS-DR 8 “main” galax¬ 
ies. 

Denoting the density fields of an LSS tracer with fgix^f), 
and that of the gamma ray emitter with /^(x, r), where r indi¬ 
cates the position in comoving coordinates and x labels time 
(given the one-to-one correspondence between time and dis¬ 
tance), the cross-power spectrum is defined as: 

(/yCu k')) = {2TTf5\k + k')P^g(k,z,z '), (2) 


where / is the Fourier transform of /(xfe)u)/(/(xfe)))> (•) 
indicates the average over the survey volume and the explicit 
dependence on z and zl highlights the possibility that the two 
populations under study ( 7 -ray emitters and extragalactic LSS 
tracers) are located at two different redshifts. From the Lim¬ 
ber approximation one gets 5{z-z'), so, in practice, only 
Pjg{k,z) is used. The modeling of the various power spec¬ 
tra used in our analysis is derived in Appendix C. Objects in 
the catalogs are described in terms of their halo occupation 
distribution (HOD), which is discussed in Appendix B. 

The window function Wg{z) appearing in Eq. (1) weights 
the contribution of objects at different redshifts to the cross¬ 
correlation signal. In the case of LSS tracers it coincides with 
the redshift distribution of the objects, dNg/dz. More pre¬ 
cisely, Wg(z) = H(z)/cdNg/dz such that f AxWgix) = 1 for a 
redshift distribution dNg/dz normalized to unity. The expres¬ 
sions of dNg/dz for the different types of LSS tracers that we 
consider here are the same as in Xia et al. (2015) (see also 
Appendix A. 3). 

For a 7 -ray emitter the window function W^ix) can be de¬ 
fined in term of the 7 -ray intensity integrated along the line- 
of-sight, If^in) as function of the direction in the sky n, which 
can be written as: 

so that = Jdx fLyCx)- We will use a coordinate system 
centered on the observer so that r = xn. The expression of the 
density fields and window functions for the different 
classes of 7 -ray sources are provided in Appendix A. 

In the Appendices we also define our reference models 
for the astrophysical and DM 7 -rays emitters, and for their 
cross-correlations with the LSS tracers. The mean intensity 

= {If_,) as function of energy of the different 7 -ray emitters 
for our reference models is shown in the left panel of Fig. 1 . 
The various curves in color indicate the contribution of each 
component, as indicated by the labels, while the black line 
indicates the sum of all astrophysical contributions. The pre¬ 
dicted total energy spectrum matches the recent Fermi-LAT 
measurements (Ackermann et al. 2015a) (solid dots with Icr 
error bars). Similarly, as shown in the right panel of Fig. 1 , we 
also verified that our reference model matches the observed 
angular power spectrum of the diffuse extragalactic gamma 
ray background measured in the 1-2 GeV energy band by the 
Fermi-LAT (grey strip) Ackermann et al. (2012). The differ¬ 
ent curves in color show the predicted angular power spec¬ 
tra of the various emitters that contribute to the total angular 
spectrum (solid black line). The model angular power spectra 
for the various gamma ray emitters have been derived by us¬ 
ing in Eq. (1) the power spectrum of the source P^^ instead 
of the cross-spectrum P^g, and W^(x) instead of the product 
WgW^. With respect to the more accurate procedure used in 
Di Mauro et al. (2014b), here we use the simplifying assump¬ 
tion that all the sources of a given population have the same 
photon spectral index (see Appendix A). 

In the next Section we will fit the theoretical predictions 
to the measured cross-correlations and will estimate the free 
parameters of the models. Eor the astrophysical components, 
the results will be in the form of deviations from reference 
models, which we adopt from the literature as updated bench¬ 
marks. We will therefore allow variations only in their nor¬ 
malization, plus a correction term as specified in the next Sec¬ 
tion. 
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Fig. 1.— Average gamma-rays intensity /-y as a function of photon energy (left) and auto-correlation APS in the (1-2) GeV energy band (right) for 

the benchmark 7 -ray models considered in this work. The black lines denote the total contribution arising from astrophysical sources (i.e. the sum of BL-Lac, 
inAGN, FSRQ and SFG emission). Fermi-LAT data are shown as black points (left, from Ackermann et al. (2015a), adding in quadrature systematic and statistical 
uncertainties) and a shaded region (right, from Ackermann et al. (2012), including the measurements with and without foreground cleaning). 


A technical remark to take into account when comparing 
the model with the data is that the experimental CAPS de¬ 
termined from the data are not deconvolved from the effect 
of the point spread function (PSF) of the instrument and the 
effect of map pixelization. To account for these effects we 
thus convolve our model prediction in Eq. (1) with the PSF 
and pixelization using the same procedure described in Xia 
et al. (2015). Formally, this is implemented defining the new 
quantity directly comparable with the data as = Vk/ 
where the effective beam window function parameterizes 
the PSF and pixelization effects (see Xia et al. (2015) for more 
details). 

Finally, in the following, we perform our analyses in terms 
of the cross correlation function CCF^'^^\6) rather than the 
cross-angular power spectra C]^. To obtain the CCF we per¬ 
form a Fegendre transformation on our CAPS as follows; 

p,[cos(0)] , (4) 

i 

where 0 is the angular separation in the sky and are the 
Fegendre polynomials. 

3. STATISTICAF ANAFYSIS 

In order to assess the possible presence of a DM signal in 
the IGRB, and its robustness to the presence of astrophysical 
emitters, we perform a statistical analysis fitting the observed 
cross-correlation data of Xia et al. (2015) with a combination 
of both DM and astrophysical source models. Specifically, 
we define a statistic from the data D, i.e., the observed 
CCF between the Fermi maps and the number of sources in 
catalogues (Xia et al. 2015), and M, i.e. the model CCF cal¬ 
culated for the different types of 7 -ray emitters as introduced 
in the previous Section and detailed in the Appendices. The 


is defined as: 


V =E E E -«r’w)) [c'-'l i, ( 




p=\ n=l 9j9j 


(5) 

where the index p runs over the five different catalogues of ex- 
tragalactic sources (2MASS, NVSS, SDSS-DR 6 QSO, SDSS- 
DR 8 Main Sample Galaxies and SDSS-DR 8 Fuminous Red 
Galaxies), the index n runs over three 7 -rays energy ranges 
(E > 0.5 GeV, £ > 1 GeV and £ > 10 GeV), whereas the in¬ 
dices 0 , and 9j run over 10 angular bins logarithmically spaced 
between 0 = 0 . 1 ° and 100 °. is the covariance matrix that 
quantifies the errors on the CCFs in each angular bin and the 
covariances among different bins, and A denotes the vector 
of free parameters which the CCF model M depends upon 
(specified below). Both the covariance matrix and the 
measured CCFs are taken from Xia et al. (2015). In 
Eq. (5) the total is obtained by adding up the individual 
X^ computed in three overlapping energy bands. There is, 
thus, in principle, a statistical dependence among the differ¬ 
ent energy bands that should be accounted for. Nonetheless, 
such dependence is expected to be small since photon counts 
are heavily dominated by events near the lower end of each 
energy interval because of the steep IGRB energy spectrum 
cx E~^'^ (Ackermann et al. 2015a). For this reason we will 
treat the CCFs estimated in the three energy intervals as sta¬ 
tistically independent in the x^ analysis. 

For any given catalog of FSS tracers, energy band and an¬ 
gular bin (i.e. for a given choice of p, n, and Oi) the theoretical 
CCF can be expressed as a sum of different contribu¬ 

tions; 


5 


M. 


(P,n) 

S, 


^A„c<f'")(0,)+A'^^4>(0,). 


(6) 
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The sum runs over the five different 7 -ray emitters: BL 
Lacs, FSRQs, SFGs, mAGNs and DM. The terms 
denote the benchmark theoretical model CCFs described in 
the Appendix and Aq is a free normalization parameter that 
quantifies the individual contribution to the observed cross¬ 
correlation. Values Aq, = 1 thus denote models equal to our 
benchmarks, while values Aq, ^ 1 would correspond to devia¬ 
tions from the benchmarks. 

Besides the normalization of each component, we have in¬ 
troduced in the fit also a further free parameter, Ai/,, which 
we dub 1-halo correction-term. This term is introduced as a 
correction for possible inaccuracies in the modeling of the 1 - 
halo contribution of the power spectrum (and thus mainly to 
the small scale cross-correlation signal) of the 7 -ray sources 
(hence its name), as discussed after Eq. (C 6 ) in Appendix C. 
For simplicity we model it as a constant term added in the 
CAPS, which is a good approximation of the 1-halo term it¬ 
self for astrophysical components, except at very high multi¬ 
poles I > 1000 , which we do not considered in our analysis. 
In real space, and taking into account the modulation intro¬ 
duced by the PSF of the instrument, the 1-halo correction- 
term explicitly reads: 

I 

where wf" is the (energy dependent) window function of 
the PSF, also introduced in the previous section. For def¬ 
initeness we have assumed that has the same energy 

dependence as the IGRB spectrum c>c . With this as¬ 
sumption we can separate the energy dependence from that 
on the source catalog where /„ = 2.46, 1,0.05 

for £„> 0.5, 1, 10 GeV (we take E = I GeV as the normal¬ 
ization energy). In this way, combining Eq.s ( 6 ) and (7) we 
have: c^"fl = f„ Notice that, in principle, each 

candidate 7 -ray emitter has its own 1-halo correction-term 
for each catalog and energy band. However they are degen¬ 
erate and only their sum can be constrained. We have thus 
grouped them together so that in Eq. ( 6 ) the reported 1- 
halo correction-term actually represents the sum of the 1-halo 
correction-terms from all the components, for a given catalog 
and energy band. 

All five values are treated as free parameters in the anal¬ 
ysis. Notice that they can be either positive or negative since 
we intend them as possible correction to our benchmarks, and 
the natural expectation would thus be A^^= 0 if the bench¬ 
marks are correct. The whole set of and A^^ coefficients 
(plus the additional parameter represented by the DM mass, 
upon which the DM signal depends) defines the parameter 
vector A of Eq. (5), which represents the full set of parame¬ 
ters over which our analysis is performed. 

Eor what concerns the particle DM contribution, we con¬ 
sider both the case where 7 -rays are produced through DM 
particle annihilation and the case of DM decay. The DM mass 
is varied from 10 GeV to 5 TeV and we will show the results 
for a DM which dominantly annihilates/decays into one of the 
following 7 -rays production channels: hh, fAfT, and 
W'^W~. Eor annihilating DM, the signal strongly depends on 
the clustering at small scales and, in particular, on the amount 
of substructures. As discussed in Appendix A, results will be 
shown for three DM substructure models: HIGH, LOW and NS. 
This will bracket the uncertainty on the reconstructed DM pa¬ 


rameters arising from DM structure modeling. The HIGH sce¬ 
nario provides a more optimistic case with the largest boosting 
factor for the 7 -ray annihilation flux. The NS scheme, where 
DM substructures are absent, provides a lower limit to to the 
annihilation signal and therefore represent the most conserva¬ 
tive scenario. Einally, the LOW scheme represents an inter¬ 
mediate case which can be (currently) considered as the most 
realistic one and that we regard as our reference model. Each 
one of these three scenarios predicts a different CCE, c^’"^( 6 >,) 
with a = DM in Eq. (6). Since the intensity of the DM sig¬ 
nal is proportional to the DM annihilation cross section {aav) 
(or DM decay rate, F^/), we normalize our calculations to the 
reference values {aav)^ = 3 • 10“^®cm^s“', i.e. the so-called 
“thermal value" which correspond to a DM particle thermally 
produced in the early Universe which, alone, would account 
for the observed DM relic abundance. For decaying DM we 
normalize the models to a decay rate of F^/ 0 = 1 -67 • 10“^^s“', 
which is the decay rate which would produce a DM signal 
equal to the one of an annihilating DM with thermal cross sec¬ 
tion in the LOW substructure scheme (for DM masses around 
100 GeV). The parameter Aq, for a = DM can be thus seen 
as the annihilation or decay rate in units of {<Jav)Q or F^o, 
respectively. 

In summary, the global fit will be performed in a 11- 
dimensional parameter space, with the parameter vector given 

by A = (Adm, fnDM,ABLLac,A„^cN,AsFC,AfSRQ,A^i^ ’ ' ’ ' 
All the parameters in the fit are linear, except moM which 
enters non-linearly in the fit through CoMiOi)- Beside the 
above fit, we will also consider different configurations, 
namely different parameter vectors A, to cross-check the 
robustness of the results. In particular, we will consider the 
case where the Aj^ are set equal to zero. These additional 
analyses will be described in more details in the next Section. 

In order to efficiently scan the multi-parameter space we 
adopt the Markov Chain Monte Carlo (MCMC) strategy pub¬ 
licly available in the cosmomc package (Lewis & Bridle 
2002). We will use linear priors limited to positive val¬ 
ues for the normalization of the astrophysical components 
ABLLacAmACNAsFcAFSRQ, although we will also check log 
priors. For the A\h parameters we allow for linear priors with 
negative values since the 1-halo correction-term can either 
correct for over-estimation or under-estimation of the small- 
scale cross-correlation. Finally we will use a logarithmic prior 
for A DM and m^M since, theoretically, the possible values of 
the DM mass and signal normalization can span several orders 
of magnitude. 

Notice that in our analysis we consider only the cross¬ 
correlation signal and ignore the intensity and the auto¬ 
correlation of the IGRB. These additional observational in¬ 
puts will be used to perform an independent a posteriori check 
on the validity of our results. While the total intensity 7^ and 

the 7 -rays autocorrelation calculated from the derived 
best-fit configurations must not exceed the measured values 
(this will be a sanity check), if they fall short of account¬ 
ing for the data this might indicate either that the measured 
IGRB contains an unaccounted contribution which does not 
correlate with extragalactic tracers (possibly of Galactic ori¬ 
gin) or that the modeling of the known components is imper¬ 
fect/incomplete. We will discuss more in detail these aspects 
in section 4.1 and in the conclusions. 

4. RESULTS AND DISCUSSION 



Cuoco et al. 



0 0.2 0.4 0.6 




0 0.2 0.4 0.6 

A 

BLLacs 


FSRQs 


Hi "1 


6 0 0.2 0.4 0.6 0 0.5 1 

A A 

SFGs mAGNs 



logjoC^DM^ 


Fig. 2.— Triangle plot of the parameters posterior distributions, for our reference fit setup, for the LOW DM substructure scheme and for a DM particle 
annihilating into bb. The darker (innermost) and the lighter (outermost) areas denote the Icr and 2a credible regions, respectively, for each combination of 
parameters considered in the analysis. The plots along the diagonal show the marginalized one-dimensional posterior distributions for each parameter. Notice 
that, for clarity, only the 5x5 sub-triangle plot with the Ac, parameters is shown, instead of the full 11x11 full triangle plot, which includes also the dark matter 
mass nioM and the 1-halo correction A® amplitudes. 


The triangle plot shown in Fig. 2 summarizes the results of 
our analysis for a benchmark annihilating DM case with bb 
hnal state and LOW substructure scheme. The plot shows the 
posterior marginal distributions of the normalization param¬ 
eters Aq. The two-dimensional plots refer to the Icr and 2a 
credible regions for each pair of parameters, while the diago¬ 
nal shows the one-dimensional posterior distribution for each 
parameter. The parameter space is eleven-dimensional, but 
for clarity we show only a part of the full triangle plot (with¬ 
out including here the parameters Aj^ and moM)- The two- 
dimensional posterior of {Aqm, moM) is shown separately in 
the left panel of Fig. 3. The posterior probability for nioM 
is instead displayed in the right panel of Fig. 3, while the 
one-dimensional posteriors for the Aj^ parameters are shown 
in Fig. 4. Finally, Fig. 5 shows the best-ht results compared 
with the measured cross correlation functions. 

A noticeable result from Fig. 2 is the fact that all the A^ 
posteriors seem to peak at A^ = 0 or close to it, except (but 
with a low signihcance) for the DM and mAGN constribu- 
tions. This does not necessarily imply that the best fit is found 
when the contributions from all components is zero. Instead, 
it is an indication that strong degeneracies are present. A two- 
dimensional analogy is given by a case with only two param¬ 


eters related by a simple relation Ai -I-A 2 = const. While the 
degeneracy would be clearly seen in the two-dimensional pos¬ 
terior, both the one-dimensional A 1 and A 2 posteriors peak at 
zero, although Al and A 2 are never both zero at the same time. 
This is precisely the results we hnd here, although the high di¬ 
mensionality of our parameter space prevents us from clearly 
trace the parameter degeneracy even in the two-dimensional 
posteriors plots. 

The degeneracy of the different astrophysical components 
can be traced to the behavior of their respective window func¬ 
tions W(z), which possess a relatively similar evolution as a 
function of redshift, and to a similar behavior of their cross¬ 
correlation 3D power spectra. As can be seen in Fig. 13 in 
Appendix A, apart from the DM case for which the 7 -rays 
emission is concentrated at low redshift with a fast decrease 
for increasing distances, astrophysical sources possess a rela¬ 
tively broad kernel. Fig. 14 and Fig. 16 instead show some 
examples of 3D cross-spectra between LSS tracers and the 
various astrophysical sources considered here or DM; we no¬ 
tice that, for a given LSS tracer (e.g. 2MASS in Fig. 14), 
the behaviors are quite similar for all astrophysical sources 
(while, instead, differences can be appreciated for the DM 
case). These facts, together with the relatively large error 
bars makes astrophysical-component separation currently dif- 
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Fig. 3.— Left: la and 2a allowed credible regions for the annihilation rate {a^v) versus the DM mass moM in the NVSS-10 = 0 (blue) and NVSS-10 
^® (tsd) fit setups. A DM particle annihilating into bb and the LOW substructure scheme are assumed. The lower Icr and 2a contours of both cases extend 
down to (ctov) = 0 (providing therefore only upper limits on the annihilation rate). Right: Marginalized ID posterior probability for the DM mass, for the same 
DM annihilation channel (bb) and substructure scheme (LOW) as in the left panel. The four lines refer to the four different fit setups described in the text, as 
labeled. 


ficult. On the other hand, given the somewhat different 3D 
cross-spectra, perspective to separate the DM component are, 
perhaps, brighter. 

An exception in this line of reasoning are mAGNs and SFGs 
which exhibit a significant degree of degeneracy with DM in 
the 2D posterior contours. The main features of the DM signal 
is that it peaks at low redshift and that is mostly contributed 
by massive halos. To mimic such a signal an astrophysical 
source must then preferentially be hosted in large halos at 
low z. Both SFGs and mAGNs meet the redshift require¬ 
ment while the blazars do not, since their window pe^s at 
higher z . However, only mAGNs are believed to be hosted in 
large halos, while SFGs typically populate galaxy-size halos. 
Objects in large halos at low redshifts are expected to have 
a large bias and, more importantly, their correlation proper¬ 
ties at the Mpc scale is dominated by a large one-halo term. 
This introduces a characteristic feature in the cross-PS that 
differentiates mAGNs from SFGs, making their contribution 
more similar to the DM one at ~ Mpc scales (see left panel 
of Fig. 15). At the lowest redshift considered (namely, in the 
cross correlation with 2MASS), the Mpc scale corresponds to 
a sub-degree scale in the CCF. Nonetheless, given the present 
still large error bars, the above feature is only weakly con¬ 
strained and thus a further degeneracy of both components 
with SFGs still remains on top of the mAGN-DM main de¬ 
generacy. Further investigation of this issue is reported later 
below. Instead, further differences between the mAGNs and 
and the DM cases are expected at smaller angles which, un¬ 
fortunately, cannot be investigated given the size of the Fermi- 
LAT PSF. 

Difficulties in modeling the 1-halo term in the HOD frame¬ 
work described in Appendix B propagates into uncertainties 
in predicting the cross-power at small-angles. To account for 
this potential source of systematic errors we introduced in 
Eq. (7) the 1-halo correction-terms The ID marginal¬ 

ized posteriors of the associated extra five parameters are 
shown in Fig. 4 as black solid curves. The various datasets 
are consistent with the case A*;, = 0 with different confidence 
levels, except NVSS. In this case we find a strong and statis¬ 
tically significant deviation from zero. This can also be ap¬ 


preciated in the fit to the observed CCF in Fig. 5 where the 
presence of a prominent 1 -halo correction-term is required to 
fit the data at small angles. There is a likely explanation for 
this additional contribution: the presence in the NVSS cat¬ 
alog of 7 -ray point sources (i.e., AGN) that are just below 
Fermi detection threshold. These sources would add their 
auto-correlation signal at zero-lad that, because of the PSF, 
spreads out to ~ 1 deg scale. This effect requires some fine 
tuning of the parameters defining the 1-halo term in Eq. (C 6 ) 
which the benchmark model fails to catch, thus requiring a 
large correction term. The effect is also discussed in Xia et al. 
(2015) to which we refer the reader for further discussion. 
The relevance of this term in the fit to NVSS data is expected 
to affect our constraints of the DM properties. To investigate 
this issue we use three further fitting procedures in addition 
to the one adopted so far. The four fitting procedures are as 
follows: 

• NVSS-10, A*^ y 0. All the 10 NVSS data points are fit¬ 
ted and the 1-halo-correction terms are free parameter 
of the fit. This is the standard fitting procedure used to 
obtain the results shown in Eig. 2. 

• NVSS-10, A^,, = 0. All the 10 NVSS data points are 
fitted and all the 1 -halo-correction terms are set equal 
to zero. 

• NVSS- 6 , A\^ y 0. The first 4 NVSS data points at small 
angles are excluded from the fit. The 1 -halo-correction 
terms are used as free parameters in the fit. 

• NVSS- 6 , Aj^ = 0. The first 4 NVSS data points are 
excluded from the fit. All 1 -halo-correction terms are 
set equal to zero. 

Eig. 4 shows the Aj^ posteriors for the NVSS-10, Aj^ 4 0 
(black solid curves) and the NVSS- 6 , Aj^, ^ 0 (red dashed 
curves) cases. It can be seen that there are no significant 
differences between the two fitting schemes, except for the 
NVSS case in which Aj^ becomes obviously unconstrained 
when the first four data points, where the fit is guaranteed by 
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Fig. 4.— Marginalized ID posterior probabilities for the I-halo correction terms A , for the NVSS-10, Aj^ yO (black-solid) and NVSS-6 , A[,j yO (red-dashed) 
fits. They are in units of 10^'^cm^^s^*. 


the 1-halo-correction term, are ignored. Fig. 6 quantifies the 
impact of the four fitting schemes on the posterior probabili¬ 
ties of all Aq, parameters. The plots show that the fitting pro¬ 
cedure does have an impact on some Aq, parameter. In particu¬ 
lar the results obtained with the NVSS-10, Aj;, = 0 fit deviates 
from the others in most of the cases. However, this is also the 
scheme that provides the worst fit to the various datasets, as 
illustrated by the comparatively larger values listed in Ta¬ 
ble 1, so that the results from this case are likely somewhat bi¬ 
ased. Instead, all the three remaining schemes provide reason¬ 
ably good fits to the different datasets. The NVSS-6, Aj^ = 0 
case provides a slightly worse fit to the data, particularly to 
the LRG sample, than the other schemes. It is interesting to 
notice that this fitting scheme favours a non-zero mAGN com¬ 
ponent (see the A^agn panel in Fig. 4) which is absorbed by 
the 1 -halo-correction term when a fitting scheme with A yo 
is adopted. This indicates that a degeneracy between the Aj^ 
and the A^acn parameter is present. We notice that in all three 
cases the is lower than total number of degrees of freedom. 
This is partly due to some unaccounted correlation between 
the three energy bins and between the different catalogues, 
and partly to the fact that the error bars are probably slightly 
over-estimated. Indeed, it is known that the algorithm imple¬ 
mented in the PolSpice software which is used in Xia et al. 
(2015) is not a minimum variance estimator of the error bars, 
i.e., it does not provide the smallest error possible (Efstathiou 


2004). 

Let us now discuss in more details the implications for the 
DM component. From the CCF plot in Fig. 5 we see that DM 
provides a significant contribution to the fit. Yet, the posterior 
probability of Aom in the bottom right panel of Fig. 2 does not 
provide a clear indication for a DM component. As discussed 
above, this is an indication that a DM signal may indeed be 
there but is degenerate with some other component, in par¬ 
ticular the mAGN one. In practice, the present datasets and 
our cross-correlation analysis cannot distinguish between the 
case of a large DM contribution with sub-dominant mAGN 
signal and that of a mAGN signal that dominates over the DM 
contribution. In fact, the situation is further complicated by 
the aforementioned degeneracy between mAGN and the 1- 
halo-correction terms. Before discussing the degeneracy issue 
more in detail, it is worth pointing out that i) our results are ro¬ 
bust to the choice of the fitting strategy, as shown in the right 
panel of Fig. 3 and in the bottom panel of Fig. 6, and that ii) 
despite the DM vs. mAGN degeneracy we are able to set con¬ 
straints on the annihilation cross-section, as shown in the left 
panel of Fig. 3, able to exclude the thermal value at 2a for DM 
masses up few tens of GeV (in the LOW substructure scheme) 
that, again, are robust against the adopted fitting scheme. We 
also verified the robustness of the DM constraints with re¬ 
spect to the choice of the priors for the astrophysical compo¬ 
nents. Specifically, we considered the case of log-flat priors 
instead of a linear-flat ones, and we found that the posteriors 
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Fig. 5.— Measured cross correlation function (CCF) (Xia et al. 2015) for E > \ GeV, as a function of the angular separation 9 in the sky, compared to the best 
fit models of this analysis. The contribution to the CCF from the different astrophysical 7 -rays emitters (BL Lac, mAGN, SFG, FSRQ) are shown by dashed 
colored lines, while their sum (“Astro Total") and the DM contribution are indicated by solid green and red lines, respectively. The 1-halo correction term is 
shown as a solid blue line. The total contribution to the CCF is given by the black solid line. The analogous plots for E > 0.5 GeV and £ > 10 GeV are shown n 
Appendix D. 


TABLE 1 

Best-fit x^f foR the four analysis described in the text, broken down into the contributions from the three energy bands (£ 05 , Ei 
AND £10 stand for £ > 0.5, 1, 10 GeV, RESPECTIVELY) AND THE FIVE CATALOGS USED. THE NUMBER OF DEGREES OF FREEDOM, VdoF, IS 
EXPRESSED AS THE TOTAL NUMBER OF DATA POINTS MINUS THE NUMBER OF FREE PARAMETERS IN THE FIT. 


Xbf 

2MASS 

SDSS-MG 

SDSS-LRG 

SDSS-QSO 

NVSS 

TOTAL 


Eos 

El 

Eio 

Eos 

El 

Eio 

Eos 

El 

Eio 

Eos 

El 

Eio 

Eos 

El 

Eio 

Eos 

El 

Eio 

AllE 

Adof 

NVSS- 10 A[^yo 

6.5 

8.5 

2.5 

4.0 

2.5 

6.3 

2.4 

2.1 

3.0 

16.8 

4.2 

6.9 

3.8 

3.7 

6.6 

33.5 

21.1 

25.3 

19.9 

150-11 

NVSS-10A[^,=0 

6.4 

12.5 

2.7 

13.5 

6.4 

8.9 

10.1 

9.5 

4.0 

13.9 

3.8 

4.9 

68.1 

84.6 

56.1 

112.1 

116.9 

76.6 

305.6 

150-6 

NVSS -6 A^^yO 

6.4 

8.8 

2.3 

3.3 

2.4 

6.8 

2.3 

2.1 

2.9 

17.4 

4.4 

7.1 

1.5 

2.1 

2.6 

31.0 

19.8 

21.7 

72.5 

138-11 

NVSS -6 A[^=0 

6.2 

11.3 

2.3 

4.8 

2.6 

6.8 

6.4 

6.3 

2.9 

19.0 

4.7 

6.2 

1.5 

2.0 

2.5 

38.0 

27.0 

20.8 

85.8 

138-6 


TABLE 2 

Best fit x^ foR the four NVSS- 6 , A„,ag]v = 0 setup, broken down into the contributions from the three energy bands (£05, £1 and £10 

STAND FOR £ > 0.5, 1, 10 GeV, RESPECTIVELY) AND THE FIVE CATALOGS USED. THE NUMBER OF DEGREES OF FREEDOM, VdoF, IS EXPRESSED AS 
THE TOTAL NUMBER OF DATA POINTS MINUS THE NUMBER OF FREE PARAMETERS IN THE FIT. 


Xbf 

2MASS 

SDSS-MG 

SDSS-LRG 

SDSS-QSO 

NVSS 

TOTAL 


Eos 

El 

Eio 

Eos 

El 

Eio 

Eos 

El 

Eio 

Eos 

El 

Eio 

Eos 

El 

Eio 

Eos 

El 

Eio 

AllE 

AboF 


7.0 

8.0 

2.3 

3.1 

2.4 

6.3 

2.2 

2.0 

3.3 

17.5 

4.3 

7.1 

1.4 

2.0 

2.6 

31.3 

18.7 

21.6 

71.6 

138-10 

A\i=0 Adm^Q 

6.0 

11.3 

2.2 

4.0 

2.7 

6.6 

6.5 

6.4 

2.8 

22.2 

5.8 

6.9 

1.5 

2.0 

2.5 

40.2 

28.3 

21.0 

89.5 

138-5 

A\^i0AoM=0 

6.9 

10.7 

3.8 

4.7 

2.2 

5.8 

2.2 

1.9 

3.4 

16.8 

4.3 

6.9 

1.5 

2.0 

2.7 

32.1 

21.1 

22.7 

75.9 

138-8 

A\i,=0Aom=0 

6.1 

14.9 

4.5 

6.6 

2.6 

6.2 

7.4 

6.3 

2.8 

19.5 

5.4 

6.8 

1.5 

2.0 

2.6 

41.1 

31.2 

23.0 

95.3 

138-3 


of the DM parameters are unaffected. Some small variations 
are present in the constraints of the astrophysical parameters, 
which is expected since at the moment the significance of the 
measurement is still not very high, and in this regime some 
prior dependence is typically still present. 

One thing to notice about the DM vs. mAGN degeneracy 
is that few mAGNs have been detected in 7 -rays so far. As 
a consequence their model contribution to the IGRB and its 
anisotropies is rather uncertain. One key quantity is the re¬ 


lation between the 7 -ray luminosity of these objects, C and 
the mass of their host halo M. Varying this relation within 
its uncertainty range, which is rather large (see e.g. Camera 
et al. (2014)), modifies the predicted cross-correlation signal. 
Fig. 15 illustrates this point. In the right panel we show the 
cross-power spectrum mAGN-2MASS galaxies (solid line) 
and how it changes when the M(C) relation is varied within its 
uncertainty band (dashed curves). Considering halo masses 
in the lower bound of the uncertainty strip significantly de- 
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Fig. 6.— Marginalized ID posterior probabilities for the Aq, terms. 

creases the amplitude of the 1 -halo term and reduces the am¬ 
plitude of the PS on Mpc scales. As a result the PS contributed 
by mAGNs will be very similar to that contributed by SFG, as 
shown in the left panel of Fig. 15, and since their window 
functions are also very similar (see Fig. 13), their contribu¬ 
tions to the cross-power become fully degenerate. 

We are therefore entitled to consider a scenario in which, 
due to this degeneracy, we set the mAGN contribution equal 
to zero and assume that it is absorbed by the SFG one. 
To explore this situation we consider four additional fitting 
schemes. In all of them we ignore the first data-points of 
the NVSS dataset (i.e. we use the NVSS -6 scheme) and set 
AmAGN = 0. The four schemes are obtained from all possible 
combination of and Aj^, that are either set equal to zero or 

let free to vary. The four combinations are explicitly shown in 
the first column of Table 2 in which we summarize the results 
of the -yA analysis. 

The inclusion of the 1-halo-correction terms improves the 
ht appreciably, although the improvement is mainly driven 
by the LRG and QSO datasets. The inclusion of DM with 
two extra-parameters also improves the ht decreasing the 
best ht by 4.3 and 5.8 for the hts with and without 1- 
halo-correction terms, respectively, with improvement mainly 
coming from a better ht to the 2MASS data. No scheme pro¬ 
vides a good ht to the CCF with the QSO for E >500 MeV. 
This is possibly an indication of an imperfect modeling of the 
energy spectrum in the QSO correlation. 

Fig. 7 shows the triangle plot for the case 7 O 


and Aj^ =0. When the mAGN contribution is suppressed 
(AmAGN =0) a non-vanishing DM component provides quite 
a good ht, with about a 2a deviation from zero. Further¬ 
more, the best ht values in Table 2 are very similar to those 
in Table 1 in which the mAGN component was included in 
the model (71.6 vs 72.5 and 89.5 vs 85.8 for the case with 
and without 1 -halo-correction terms, respectively), a fact that 
corroborates the evidence for a degeneracy between DM and 
mAGNs. Fig. 8 illustrates the robustness of these results to 
the inclusion of the 1 -halo-correction term, A\i^ 7 ©, for two 
different hnal annihilation states; bb (left set of plots) and and 
(right plots). In all cases the best hts preference is for 
the presence of a DM component, even when the extra degree 
of freedom Aj^ is included. Furthermore, these plots reveal 
a degeneracy between Aom and moM which is to be expected 
since, as can be seen in Eq. (Al), the WIMP signal approx¬ 
imately scales with Adm/i'^dm- Indeed, Wgi contains a fac¬ 
tor (cTflV) jinDM^ plus an integral over the energy which intro¬ 
duces a contribution roughly proportional to nioM (being the 
spectrum integrated up its endpoint, which is iuom)- Inspec¬ 
tion of Fig. 7 also shows more clearly the degeneracy between 
DM and the SFG components, left from the mAGN-DM-SFG 
degeneracy after removing the mAGN component. One con¬ 
sequence of this is that performing a ht excluding either the 
SFG or the mAGN component would enhance the strength of 
the DM signal, without affecting appreciably the value of the 
best ht 
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Fig. 7.— Triangle plot for the NVSS-6, = 0, = 0 fit. The results refers to the bb annihilation channels and the LOW DM substructure scheme. 
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refers to the bb annihilation channels and the LOW DM substructure scheme. Right: The same as in the left panel, but for the annihilation channel. 


Fig. 8 .— Left: Detail from Fig. 7 showing the DM parameters only. Furthermore, the ID posterior panels show both the Aj^ = 0 and ^0 cases. The plot 


In conclusion, our analysis indicates that a significant DM 
contribution to cross-correlation is entirely plausible. How¬ 
ever, the degeneracy with other astrophysical sources, namely 
mAGNs and SFGs, largely originating from the current obser¬ 
vational uncertainties, prevents us from drawing a dehnitive 
conclusion. Future analyses with increased 7 -ray statistics 
and improved angular resolution in which the cross correla¬ 


tion is extended to other catalogues of extragalactic objects 
will help to break the present degeneracies and to pinpoint the 
correct scenario. 

As a hnal remark, we also mention that, as a cross-check, 
we have performed the analysis employing the astrophysical 
models adopted in Xia et al. (2015). The constraints on the 
7 -ray astrophysical contributions are different, which is ex- 
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Fig. 9.— Comparison of the posteriors distributions for the DM parameters (mass widm (left) and annihilation rate in terms of the thermal one Adm (right)) for 
the NVSS-6, = 0, A,„yigjv = 0 fit and for the 3 different SFG models described in the text. 


pected given the different modeling. Regarding DM, using 
the NVSS- 6 , A^agn = 0, = 0 fit configuration, which is 

the closest to the one used in Xia et al. (2015), we compare in 
Fig. 9 the DM posteriors from 3 different fits using 3 different 
SFG models: the one adopted in this work (black solid curve), 
the SFGl model from Xia et al. (2015) (red dashed curve) 
and a modified version of SFGl (blue dot-dashed curve) with 
redshift-dependent bias equal to the the bias of the present 
model (while the original SFGl model has bias equal to 1 
for all z). The SFG2 model of Xia et al. (2015) is very sim¬ 
ilar to present SFG model and is not considered. The plot 
indeed shows that the DM results are not significantly depen¬ 
dent from the SFG model adopted. 

With no unambiguous indication for a DM components we 
can nevertheless set constraint on the properties of the DM 
candidates. To this purpose we perform, for any given mass of 
the DM particle candidate, an individual 10 parameter fit and 
set the 95% bound on (uav) from the posterior distribution. 
The results for the annihilating DM are summarised in Fig. 10. 
In the left panel we focus on the bb annihilation channel. The 
solid line with different colours refer to constraints obtained 
from each of the three energy band separately {E > 0.5,1,10 
GeV) as well as the ones obtained by their combination (red 
line). All theee results refer to the LOW substructures model 
and are obtained with the reference NVSS-10 A*,, =^0 fitting 
scheme. The upper and lower dot-dashed curves show how 
the bounds change in the the NS and HIGH substructure model, 
respectively. 

In the right panel we compare the results of different final 
annihilation channels and W^W~) to the original 

bb case (red curve) shown in the left plot, for the LOW sub¬ 
structures scenario and combining all energy bands. All the 
results refer to the benchmark NVSS-10 Aj^ case, but the 
other fitting schemes provide nearly indistinguishable con¬ 
straints. The black curve is taken from Regis et al. (2015) and 
refers to the case in which we assumed that all the 2MASS 7 - 
ray correlation is produced by DM, with no astrophysical con¬ 
tribution. As expected, including the astrophysical sources 
makes the constraints stronger, of about a factor of 4. The gain 
is significant and will further improve once the DM-mAGN- 


SFG degeneracies discussed above will be removed. 

As expected, uncertainties on the bounds driven by the sub¬ 
structure model are significant. The left panel of Fig. 10 
shows that assuming the HIGH model would strengthen the 
constraints on the cross section by about one order of magni¬ 
tude, whereas in the NS scenario, the bounds would weaken 
by about a factor of 5. This implies that the thermal annihila¬ 
tion rate {aav) = 3 • 10“^®cm^s“* is excluded at the 95 % level 
up to masses of 6 , 25, 250 GeV in the NS, LOW and HIGH 
scenarios, respectively. 

In Fig. 11 we instead show the 95% lower bounds on the 
lifetime of a decaying DM particle, for various decay final 
states. Bounds on DM decay, being proportional to the DM 
density (and not DM density squared, as instead the annihila¬ 
tion signal) depend on the total DM mass in structures and are 
not affected by the different substructure modeling. As for 
the annihilation case, including the astrophysical sources in 
the analysis improves the constraints, again by about a factor 
of 4, with respect to those obtained by ignoring the astrophys¬ 
ical components (Regis et al. 2015). 

Finally, to test the robustness of our DM constraints we 
have repeated the analysis using the same astrophysical mod¬ 
els used in Xia et al. (2015) and we found that they are very 
similar to the ones obtained in the present analysis. 

4.1. Self consistency tests: mean intensity and 
auto-correlation of the IGRB 

As anticipated in Section 3 instead of including the mean 
IRGB intensity and its auto-correlation in the fit, we use 
these additional observational inputs a posteriori as a self- 
consistent test for our best fitting model. 

We define, the fractional mean IGRB intensity pre¬ 

dicted by the cross-correlation fit, as follows 

hoT^IGRB-^FSRQhsRQ +ABLLaJBLLac + 

AmAGNlmAGN~^^SFGlsFG~^^DMlDM ) ( 8 ) 

where are the integrated 7 -ray intensities of our refer¬ 
ence models for the five 7 -ray emitters considered here and 
shown in Fig. 1 and n= 1,2,3 identifies the energy band The 
total intensity is defined as Ijqj = where the sum 
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mDM[GeV] 7nDM[GeV] 

Fig. 10.— Left: 95% upper bounds on the DM annihilation rate (cTav) as a function of the DM mass, for the LOW substructures model and the reference 
NVSS-10 tO fit. Solid lines refer to the bb annihilation channel: the red line refers to the analysis that combines information from all the three energy bins 
under consideration {E > 0.5,1,10 GeV), while the other three lines refer to the analysis performed on a single energy bin (as stated in the figure label). The 
upper dot-dashed blue line refers to the NS substnacture model, while the lower dot-dashed black line to the HIGH substructure model. Right: in addition to the bb 
case (red line) reported in the left panel, the different lines show the upper bounds for the (blue), (green) and W^W~ (magenta) annihilation channels, 
for the LOW sub-structures model. The black line instead shows the upper bound for the bb case and LOW substructure scheme, obtained under the assumption 
that the DM contribution to the 2MASS cross-correlation is the dominant one (taken from Regis et al. (2015)). 



mDM[GeV] 

Fig. 11.— For a decaying DM, 95% lower limits on the DM lifetime r as 
a function of its mass, for different decay channels: bb (red), (blue), 

r"^T~ (green) and W'''W~ (magenta). The black line instead shows the lower 
bound for the bb case obtained under the assumption that the DM contribution 
to the 2MASS cross-correlation is the dominant one (taken from Regis et al. 
(2015)) 


runs over the five types of emitters. In our model Ijqj= 
10“®, 4 X 10“^, 1.5 X sr“* for the energy ranges 

E > 0.5, 1, 10 GeV, respectively, which are consistent with 
the measured IGRB (Ackermann et al. 2015a). We thus ex¬ 
pect that the have values close to unity to match ob¬ 

servations. Note that the parameters need not to be the 
same in each energy band since the total signal Ij^j and the 
individual contributions have different scaling in energy. 
However, the difference is not very large. 

Similarly, we dehne the IGRB auto-correlation predicted 
from the cross-correlation fit, as a fraction of the measured 


one, in terms of the parameter fcpjoRB 

CpjorfcpjGRB =^FSRQ^P,FSRQ +J^BLLacCp.BLLac + 

A^mACNCp.mAGN+A^SFGCp.SFG > (9) 

where Cp,fsrq = 1.6 x 10“'*, Cp^blluc = 7.9 x 10“**, Cp^mAGN = 
3.9 X 10“*® and Cp^sFG = 6.3 x 10“^*, all of them in units of 
(cm“^s“*sr“*)^sr. are the predicted average auto-correlation 
signals in the multipole range i = 155-504 and in the energy 
band 1-2 GeV. We have neglected the DM contribution since 
it is largely subdominant with respect to FSRQs, BLLacs and 
mAGNs (see Fig. 1). The SFG contribution, which is also 
subdominant, is considered for the sake of completeness. In 
the above equation we made the assumption that the ampli¬ 
tude of the auto-correlation signal scales with the square of 
the normalization parameters of the individual components. 
Unlike the cross-correlation case we did not include any 1- 
halo-correction term since we model astrophysical emitters 
as point sources for which no additional small-scale power 
is expected to contribute to the auto-correlation signal. Like 
the mean intensity, the value /crjorb = 1 characterizes a model 
which saturates the measured IGRB auto-correlation. 

In Fig. 12 we show the posterior probabilities for A"^^^ in 
the three energy bands considered in our analysis, and fcpjoRB ■ 
We hnd that the typical value Ajcrb is between 20% and 50% 
in the two lower energy bands (upper panels) whereas for 
£ > 10 GeV is in the broader range 10% - 80%. These results 
are robust to the details of the fitting procedure, as demon¬ 
strated by the similarity of the various curves. They imply that 
the extragalactic sources considered in our model (BL Lac, 
mAGN, SFG, FSRQ and DM) which, as we have seen, pro¬ 
vide a good match to the observed cross correlation with LSS 
tracers, also contribute to a significant fraction of the IGRB, 
although possibly not to the whole signal. This result is in¬ 
teresting but should also be taken with a grain of salt given 
the complexity of our cross-correlation model. For example 
Xia et al. (2015), using a different model for the SFG emis- 
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Fig. 12. — Marginalized ID posterior probabilities for the cumulative fractional contribution Ajqkb of all 7 -ray sources (BL Lac, mAGN, SFG, FSRQ and DM) 
to the total intensity Iiqrb measured by Fermi. A/qbb is expressed in terms of Iiqbb = 10“^, 4 X 10^^, 1.5 X 10“* cm“^ s“* sr“* for the energy bins E > 0.5, 1,10 
GeV (to account for spectral behaviour). The bottom right panel show the same information for the IGRB angular auto-correlation in the 1-2 GeV energy band. 
The various lines refer to the four fits described in the text with same color and dashing conventions of Fig. 6 . 


sion and bias, was able to account for a larger fraction of the 
IGRB, although again not 100%. It should be also noted that 
the measurement of the IGRB in Ackermann et al. (2015a) is 
affected by systematic errors induced by the imperfect model 
of the foreground Galactic emission, even if the size of this 
systematic uncertainty does not seem to be large enough to 
saturate our models to 100% of the total emission. If indeed it 
turns out that additional 7 -ray sources are required to explain 
the total intensity of the IGRB, then the results of our anal¬ 
ysis set a rather sever constraint: their correlation with LSS 
tracers must be weak. This would imply that they should be 
local, possibly of Galactic origin, like the millisecond pulsar 
or, perhaps, diffuse inverse Compton photons from cosmic- 
ray electrons scattering on the optical/infrared Galactic inter¬ 
stellar radiation held. Future analyses with newer and addi¬ 
tional datasets will help to clarify this interesting issue. 

The posterior for fcpjaRs instead consistent with unity, al¬ 
though its probability distribution actually spans several or¬ 
ders of magnitude from 10 “^ to 10 , meaning that the mea¬ 
sured auto-correlation does not provide a very stringent cross¬ 
check. 

5. SUMMARY AND CONCLUSIONS 


In this paper we have used the cross-correlations recently 
measured in Xia et al. (2015) between Fermi-LAT diffuse 7 - 
ray maps and different catalogs of LSS-tracers to investigate 
the origin of the IGRB and the nature of the various sources 
that may contribute to it, including DM annihilation or decay. 
This work extends that of Regis et al. (2015) which used only 
the 7 -ray- 2 MASS correlation and considered DM as the only 
source of the extragalactic 7 -ray signal. Our main results are 
as follows: 

• Our theoretical models provide a good ht to the cross¬ 
correlation measured in all employed catalogs of extra- 
galactic tracers, namely SDSS-DR 6 quasars, 2MASS, 
NVSS, SDSS-DR 8 LRGs and SDSS-DR 8 MG. The 
quality of the ht is quantihed by means of a 7 ^ analysis 
in which we account for covariance among the errors 
in different angular bins whereas we ignore the covari¬ 
ance among energy bins and among the different cat¬ 
alogs. The hrst approximation is justihed by the pho¬ 
ton statistics, which is dominated by low energy event, 
making each of the energy bins considered in our analy¬ 
sis effectively independent. The second approximation 
is justihed by the spatial distributions of the objects in 
the different catalogs that, with the partial exception of 
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the NVSS one, do not signihcantly overlap with each 
other. 

In our cross-correlation function (CCF) models we con¬ 
sider four different types of astrophysical sources (two 
flavours of blazars, FSRQs and BL Lacs, SFGs and 
mAGNs) and, in addition, annihilating/decaying DM. 
The rationale behind the choice of these astrophysical 
sources is that previous analyses have shown that they 
are the main contributors to the IGRB and its angular 
auto-correlation. These two observational constraints 
are not considered in our fit. Instead, we use them a 
posteriori to check the consistency of our best fitting 
models which are based solely on the measured CCF. 
We hnd that models that provide a good match to the 
cross-correlation fall short of accounting for the mean 
7 -ray intensity. The discrepancy is not large, less than 
a factor of two, especially in the high energy band, and 
could be accounted for by a combination of model un¬ 
certainty and imperfect subtraction of the Galactic fore¬ 
ground. However, it may also indicate that additional 
types of sources that do not cross-correlate with the 
LSS, like 7 -ray sources within our Galaxy, are required 
to account for the whole IGRB intensity. 

Including DM among the possible IGRB sources does 
not significantly improve the quality of the ht, and does 
not indicate a preference for a particular DM mass or 
annihilation cross-section/decay rate. We hnd that the 
reason for the low statistical signihcance on the pres¬ 
ence of a DM component does not lie in the fact that the 
ht rejects this component, while it is rather due to the 
presence of a model degeneracy with other types of as¬ 
trophysical sources, mainly mAGNs and SFGs. In other 
words, a signihcant DM contribution gives an equally 
good ht as a case with a negligible DM contribution and 
a larger mAGNs and SFGs emission. Neglecting the 
mAGN component in the ht partially breaks this degen¬ 
eracy and provides a small 2a) preference for DM. 
The best ht is found for a rather canonical WIMP DM 
candidate with itidm ^ 100 GeV that annihilates into 
bb at a rate which is of the order of the thermal value 
for the benchmark LOW DM clustering scenario consid¬ 
ered. A candidate with a slight smaller mass of about 
30 GeV that annihilates into t~ provides an equally 
good ht. 

Breaking this degeneracy is the main goal of future 
cross-correlation analyses similar to the present one. 
Fortunately, this is a realistic goal. One of the main rea¬ 
son for this degeneracy is the uncertainty on the mAGN 
and SFG luminosity in 7 -ray which, to date, has been 
directly measured for a handful of very nearby objects. 
However, their number is bound to increase thanks to 
the fact that Fermi-LAT will keep taking data in the 
next few years. In addition, the quality of the Fermi 
maps is also expected to increase both in terms of pho¬ 
ton statistics, which will allow to better sample the en¬ 
ergy behaviour and to improve the sensitivity to char¬ 
acteristic DM spectral features, and angular resolution, 
which would allow us to push the correlation analysis 
to smaller angular scales where the 1 -halo term domi¬ 
nates. 


nihilation cross-section/decay rate as function of the 
DM mass. Our derived constraints are comparable 
in strenght to most of the current indirect detection 
method that exploits the 7 -ray sky (Regis et al. 2015). 
These constraints are rather robust to the astrophysi¬ 
cal details of the models but, as expected, do depend 
on the detail of the DM substructure and small-scale 
clustering. For this reason and with the aim of bracket¬ 
ing current theoretical uncertainties, in addition to the 
LOW scenario which represents the current, somewhat 
conservative, benchmark substructure model, we have 
explored two additional, rather extreme cases: the NS 
case in which we completely ignore substructures and 
that provides extremely conservative constraints of the 
DM properties, and the HIGH scenario in which sub¬ 
structures are more numerous and have an higher den¬ 
sity concentration. In the most conservative NS sce¬ 
nario our method excludes, at a credible level larger 
than 95%, that DM particles with masses smaller than 
10 GeV annihilating entirely into bb could have a ther¬ 
mal cross section. In the optimistic scenario, the same 
statement applies to particles lighter than ~ 600 GeV. 
The bounds are a factor of ^ 4 stronger than the most 
conservative case considered in Regis et al. (2015) in 
which only DM is used in order to saturate the 2MASS 
cross-correlation. Constraints on DM decay time for 
DM decaying into bb are ~ 10^* s, roughly indepen¬ 
dently from the DM mass. 

All in all we are confident that the results obtained in our 
analysis, which are already quite remarkable considering that 
this is the first time that a genuine cross-correlation signal is 
detected in the Fermi-LAT 7 -ray maps, will soon improve sig¬ 
nihcantly. In this respect this work also represents a proof of 
concept that illustrates the potential of the cross-correlation 
analysis. We base our optimism on the fact that the new 
Pass 8 data, with improved effective area and angular reso¬ 
lution, will soon be released by the Fermi-LAT Collaboration 
and that additional catalogs of objects al relatively low red- 
shifts with wide, almost all-sky, angular coverage and well 
determined redshift distribution are already available (Bilicki 
et al. 2014) and some new ones are being compiled (Bilicki 
2015). 
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APPENDIX 

WINDOW FUNCTIONS 

In this Appendix we discuss the modeling of the window functions adopted for the calculation of the cross-correlation angular 
power spectrum C'''® of Eq. (1), which are in turn the ingredient for the determination of the cross correlation function CCF^^{9) 
defined in Eq. (4). 


Dark matter 
Annihilating dark matter 

DM annihilations in haloes and in their substructures produce 7 -ray photons. This emission traces the DM density squared 
therefore the density held responsible for the correlation signal is fs 2 {x^r) = PomCXj*")- The window function reads; 


Ws2(x)-- 


(flpMPcf {fJgV) 
Att 2mnM^ 


[ 1 +z(x)]'A 2 (x) 




J E-y> Em 


ANg 

AE^ 




(Al) 


where Udm is the cosmological abundance of DM, pc is the critical density of the Universe, otdm is the mass of the DM particle, 
and {agv) denotes the velocity-averaged annihilation rate, assumed here to be the same in all haloes. AtNg/AE^ indicates the 
number of photons produced per annihilation event, and sets the 7 -ray energy spectrum. We will consider annihilation into bb 
quarks as representative of a typical soft annihilation spectrum (with 7 -rays mostly arising from production and decay of neutral 
pions), and into pApT leptons as representative of a hard-spectrum channel (where 7 -rays mostly arising from hnal state radiation), 
with and W'^W~ hnal states as intermediate possibilities. Emin is the energy threshold of the Fermi-LAT maps considered 
in the analysis, namely: Emin = 0.5, 1,10 GeV. The factor exp{-T[x,E-y(x)]} accounts for absorption due to the extra-galactic 
background light, and we model the optical depth r as in Franceschini et al. (2008). 

A crucial quantity in Eq. (Al) is the so-called clumping factor A^(x): 


A 2 (z) ^ 


'dm) 
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'dM^(M,z)[l + hsub(M,z)] 
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d^x 


pI(x\M,z) 


Pdu 


(A2) 


The clumping factor involves the integral of the halo number density An/AM above the so-called minimal halo mass Mmin, 
multiplied by the total number of annihilations produced in the generic haloes of mass M at redshift z with density prohle 
p/,(x|M,x) and with subhalos providing a “boost” to the emission given by b^uh- We assume a reference value of IQ~^Mq for 
4^min, which coiTesponds to a typical free-streaming mass in the WIMP DM scenario. We adopt the halo mass function from 
Sheth & Tormen (1999) and we assume that the halos are characterized by the so-called Navarro-Frenk-White (NEW) universal 
density prohle (Navarro et al. 1997). The prohle is completely determined by the total mass of the halo and by its size. We 
express the latter in terms of the concentration parameter c(M,z), taken from Prada et al. (2012) (see also Sanchez-Conde & 
Prada (2014) for an analytic ht of c(M, z = 0) of Prada et al. (2012)). 

Concerning the boost provided by subhaloes hosted in the main haloes, we consider three scenarios (HIGH, LOW, and NS) as 
extreme cases bracketing the effect. In Eq. (A2) this is indeed the most uncertain quantity (Sanchez-Conde & Prada 2014; Gao 
et al. 2012; Fornasa et al. 2013; Ng et al. 2014). In the LOW scenario, the function bsub{M,z) is computed following Sanchez- 
Conde & Prada (2014) (see, in particular, their Eq. (2) assuming a subhalo mass function dn/dM^uh oc M//^/). The HIGH scenario 
stems instead from the bsub{M,z) found in Gao et al. (2012), and assuming no redshift dependence. In the NS case, we simply set 
^sub = 0 ; this can be considered as the most conservative approach. 

The blue curves in Fig. 13a show a quantity related to the window function of Eq. (Al), dehned as the kernel of the 7 -ray 
emission entering in the computation of the angular power spectrum discussed in Section 2. Specihcally, we dehne the kernels 
as: 

Kg = \/cz/Hiz)Wg(z)/xi^) , (A3) 


such that; 


AiiSi) 


d \niz)K^, (z) Kg^ (z) P^,g. (,k,z). 


(A4) 


In Fig. 13a we choose a reference particle-physics model with todm = 100 GeV, {agv) = 3 x 10“^®cm^s“' and bb annihilation 
channel. The three clustering scenarios (NS, LOW and HIGH for dotted, solid and dashed curves, respectively) share approximately 
the same redshift dependence, but they correspond to different sizes of the clumping factor and consequently of the intensity of 
the DM-induced 7 -ray flux. Note that a comparison with previous works in the literature can be non-trivial, as different groups 
employ different prescriptions for the ingredients of the DM clustering and, in particular, for the boost factor. Fig. 13 can be 
useful also as a normalization test when confronting the results presented in the rest of the paper with other works. 


Decaying dark matter 

If instead of being stable, the DM particles decay, while having a negligible self-annihilation rate, the produced 7 -rays traces 
the DM density linearly, i.e., fsix^^) = Pdm(X)*‘), The window function in this case reads; 
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Fig. 13.— Angular-power-spectrum kernels K(z) of 7 -ray emitters (left) and galaxies (right), shown as a function of the redshift z- The kernel is defined as 
Ka(z) = y/cz/H(z)Wa{z)/x(z)^ where H{z) is the Hubble parameter, Wa(z) denotes the window function of the objects of class a and x is the comoving distance. 
The 7 -rays kernels are integrated for energies above 1 GeV and refer to unresolved sources fainter than Fsens = 5 • 10“^^ photons cm“^ s“*. 

where Tj is the decay rate. The photon yield, dN^/dE^, is assumed to be the same as for annihilating DM, but with the energy 
of the process given by ,/s = todm instead of 2 otdm- In other words, dNd/dE/E/ = dtNa/diE/2E/) with the kinematic end-point 
being at The kernel in the case of decaying DM is shown as a cyan curve in Fig. 13a. In the plot we report reference 

particle-physics model with mDM = 200 GeV, Td = l/Fd = 6 x 10^^ s and decays into hh quarks. Note that for decaying DM, the 
window function does not depend on the details of the DM clustering. We notice also that DM kernels peak at low redshifts, both 
for annihilating and decaying DM, and have a relative fast decrease with distance. 


Astrophysical sources 

For astrophysical sources, we adopt as the characterizing parameter the source 7 -ray luminosity C in the energy interval (0.1 - 
100) GeV. For a power-law energy spectrum with spectral index a, the window function takes the form; 


Ws,{x) 




(A 6 ) 


where Eq = 100 MeV is just the normalization energy, and i stands for each of the 7 -rays sources adopted in our analysis: BL 
Lac, FSRQ, mAGN and SFG. The mean luminosity produced by an unresolved class of objects located at a distance x from us is 
denoted by (/s,(x)) and is given by: 

/ •^max (Fsens 5^) 

dCC-^iiCx), (A7) 


where $,(£,z) is the 7 -ray luminosity function for the source class i. The upper bound, >Cmax(Fsens,z). is the luminos¬ 
ity above which an object can be resolved, given the detector sensitivity Fjens for which we assume the value Ljens = 5 x 
10 “'° photons cm“^s“* above 1 GeV (Nolan et al. 2012 ; Acero et al. 2015 ). The precise value depends slightly on a, and on 
the catalogue of resolved point sources, although varying Fsens within these different values has only a weak impact of the win¬ 
dow function. Conversely, the minimum luminosity Cmin.i depends on the properties of the source class under investigation. The 
four populations of astrophysical 7 -ray emitters (i.e., BL Lac, FSRQ, mAGNs and SFGs) are discussed in the following. For 
each of them we describe the choice of a, and of the 7 -ray luminosity function. 


Blazars 

We consider BL Laceitae (BL Lacs) and flat-spectrum radio quasars (FSRQ) separately. The 7 -ray luminosity function of BL 
Lacs and FSRQ is taken from Ajello et al. (2014) and Ajello et al. (2012), respectively, where it is derived from a parametric 
fit of the redshift and luminosity distributions of resolved blazars in the Fermi-LAT catalogue. The lower limit of the integral 
in Eq. (A7) is set to Emin = 7 • 10"'^ergs“* (BL Lac) and Emin = 4 • 10''^ergs“* (FSRQ). For the energy spectrum, we consider a 
simple power-law with a spectral index taken from the average spectral index in Ajello et al. (2014, 2012), namely, we assume 
ccfiLLac = 2.1 and ctfsrq = 2.44. 

The kernels of unresolved blazars are shown by the solid red (BL Lac) and magenta (FSRQ) lines in Fig. 13a. Note that they 
strongly decrease at low z since Fermi-LAT has already detected a large number of the closest (brightest) emitters of these classes. 
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Misaligned AGNs 

In the case of mAGN, we follow Di Mauro et al. (2014a), which studied the correlation between the 7 -ray luminosity and the 
core radio luminosity Lr.core at 5 GHz, and derived the GLF from the radio luminosity function. We consider their best-ht C vs. 
Lr.core relation and assume an average spectral index amACN of 2.37. The solid green line in Fig. 13a indicate the contribution of 
unresolved mAGNs. 


Star-forming galaxies 

As done in Ackermann et al. (2012), we assume that the 7 -ray and infrared (IR) luminosities are correlated in the case of SFG. 
We adopt the best-ht C vs. Ljr relation from Ackermann et al. (2012) while for the IR luminosity function we adopt the one 
from Gruppioni et al. (2013) , (adding up spiral, starburst, and SF-AGN populations of their Table 8 ), as considered in Tamborra 
et al. (2014). The spectral index is taken to be osfg = 2.7 for all the 3 components although starbursts galaxies would require in 
principle a somewhat harder spectrum. Nonetheless, this component is subdominant in the total SFG contribution except for high 
energies and at high redshift (i.e., in the ranges which are less relevant for the analyses in our work). The above choice has thus 
no practical effects on our results. The kernel associated to unresolved SFGs is the solid orange line in Fig. 13a. All the different 
single peaked sub-populations provide sizable contributions and this gives raise to different peaks. 


The 7 -ray emission produced by the four extragalactic astrophysical populations described above accounts for approximately 
the whole IGRB and autocorrelation angular power spectrum (see Fig. 1). As described in the main text we however introduced 
a normalizing constant Aq, for each population to be determined by the fit. Apart from the extragalactic DM-induced emission 
described in Secs A. 1.1 and A. 1.2, there may be a contribution associated with annihilations/decays in the DM halo of the Milky 
Way. This is not included since it does not correlate with the LSS tracers. 

Galaxy catalogues 

For galaxies, we take the redshift distributions dNj/dz(x) reported in Xia et al. (2015). The associated kernels are shown in 
Fig. 13b. The 2MASS kernel peaks at low-redshift, and a comparison with the 7 -rays kernels shown in the left panel of the same 
Fig. 13 indicates that the 2MASS catalogue is the most suitable for investigating a DM signal in the cross-correlation analysis, 
followed by the SDSS Main Galaxy Sample catalogue. 


HALO OCCUPATION DISTRIBUTION OF GALAXIES 


In this work, we compute the angular cross-correlation between the unresolved 7 -ray sky and the number of galaxies in specihc 
catalogues. In order to estimate the latter from a theoretical point of view (and since we adopt the halo model description for the 
structure clustering), we need to describe how galaxies populate halos. Namely, we need to model how many galaxies of a certain 
catalogue are present in a halo of mass M and how they are spatially distributed. To this aim, we employ the halo occupation 
distribution (HOD) formalism. 

We follow the approach described in Zheng et al. (2005) (for review on HOD, see also Berlind & Weinberg (2002); Cooray & 
Sheth (2002)), where the HOD is parameterized by distinguishing the contributions of central and satellite galaxies, N = Ncen+Ns^t 
(since different formation histories typically imply different properties for galaxies residing at the centers of halos with respect 
to satellite galaxies). These can be modeled with the following functional forms; 


(Neen(M)) = 
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With this formalism, we need hve parameters for each galaxy population; Mth denotes the approximate halo mass required to 
populate the halo with the considered type of galaxies, with the transition from 0 to 1 central galaxy modeled by means of 
Eq. (Bl), and set by the width cTLogM- The satellite occupation is described by a power law (with index a and normalization 
set by the mass Mi), with an exponential cutoff Mcut at low masses. The value of the hve HOD parameters for each of the 
considered galaxy population is discussed in the following. Eor some catalogues, we will also consider similar but slightly 
different functional forms. 

We selected those galaxy samples with available HOD which more closely resemble the catalogues considered in the cross¬ 
correlation analysis of this work. We caution, however, that since the matching of the two samples is not perfect some differences 
in the associated HODs might be expected. Nevertheless, this should not affect our results in a dramatic way. 

Eqs. (Bl) and (B2) provide the number of galaxies in a halo of mass M. Concerning the spatial distribution, we treat central and 
satellite galaxies separately. The former is taken as a point-source located at the center of the halo (the point-source approximation 
is expected to break down only for £ > 10^). Satellite galaxies are instead described in an effective way with a spatial distribution 
following the host-halo prohle. In other words, we express the density field of galaxies with; 


gg(x-x'\M) = (Ncen(M)) 5^{x-x')-\- {N^^M)) ph{x-x'\M)/M. 


(B3) 


/ 


d^ 


Xggix) = (Neen(M)) + (M,t(M)) = (A(M)) . 


Note that; 


(B4) 
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A determination of the HOD for the 2MASS galaxies is not present in the literature (to the knowledge of the authors). In 
Zehavi et al. (2005), a sample of about 200,000 SDSS galaxies mostly residing in the redshift range 0.02 < z < 0.167 and with 
r-band magnitude 14.5 <r< \1.11 was analyzed. Such ranges of redshift and magnitude are analogous to the ones of the adopted 
2MASS catalogue (Jarrett et al. 2000), with the cross-identihcation of the latter with SDSS found to be successful for about 90% 
of the sources (McIntosh et al. 2006). We can thus exploit the HOD results of Zehavi et al. (2005). They considered a step 
function ({Ncen) = 0 for M < Mth and (Ncen) = 1 for M > M±) instead of Eq. (Bl), and set Mcut = 0 in Eq. (B2). The analysis 
was performed by splitting the sample in luminosity bins, but for our purposes we can consider the averaged best-ht parameters 
weighted over the number of galaxies in each bin. We found logM* = 12.1, a = 1.2, and logMi = 13.5. 

NVSS HOD 

The NVSS sub-sample considered in this work (Blake & Wall (2002)) contains sources brighter than 10 mJy at 1.4 GHz. The 
vast majority of them is associated to bright AGNs. To model the AGN HOD, we follow Chatterjee et al. (2012). Eor bright 
objects, they found logM* = 13.03, CLogM = 0.96, a = 1.17, \ogMcut = 11.5, and logMi = 13.64. 

SDSS-DR8 Main Galaxy Sample HOD 

In Ross et al. (2010), the clustering of more than three million photometrically selected SDSS galaxies was analyzed. In 
particular, the sample was defined requiring de-reddened r-band magnitudes rj < 21 and absolute magnitudes Mr < - 21 . 2 , in the 
redshift range 0.1 < z < 0.4 and masking objects with Galactic extinction > 0.2. This galaxy sample is very similar to the one 
adopted in this work for the cross-correlation analysis (Aihara et al. (2011)), except for a more limited redshift (and magnitude) 
range. Since the peak of the redshift distribution is at z ~ 0.3, such difference is not expected to play a major role. After averaging 
the best-fit values of HOD parameters over the different redshift bins considered in Ross et al. (2010) (weighted for the number 
of galaxies in each bin), we obtain logMth = 12.09, (TLogM = 0.3, a = 1.09, and logMi = 13.25. The functional form of the satellite 
HOD considered in Ross et al. (2010) is: 

(N™,) = (Nce„)x f (B5) 

instead of Eq. (B2), with (Veen) from Eq. (Bl) and (Vsat) = 0 forM < Mth. 

SDSS-DR8 Luminous Red Galaxies HOD 

A recent analysis of more than 500,00 SDSS-III CMASS galaxies (Reid et al. 2014) derived the HOD of this galaxy sample. 
The satellite HOD was modeled by means of: 

(Vsat) = (Veen) X (B6) 

with (Veen) given by Eq. (Bl) (and (Vsat) = 0 for M < Meut). The best-ht parameters in the hducial model of Reid et al. (2014) 
are logM* = 13.03, CLogM = 0.38, a = 0.76, logMi = 14.08, and Meut = 13.27. These results are found to be in agreement with 

the HOD analysis presented in White et al. (201 1), which in turn was tested to reproduce the clustering of the galaxy sample of 

Ho et al. (2012) adopted here. 


SDSS-DR6 Quasar HOD 

The modeling of the halo occupation distribution of SDSS quasars is taken from Richardson et al. (2012). The sample consists 
of 47,699 quasars in the redshift range 0.4 < z < 2.5 with median redshift of z= 1.4 and Hux limited to i < 19.1. It is very 
similar to the catalogue considered in this work for the cross-correlation analysis (Richards et al. (2009)). The best-ht parameters 
entering in Eqs. (Bl) and (B2) are not provided in Richardson et al. (2012), but we hnd that with logMth = 16.5, (TLogM = 165, 
a = 1, logMcut = 15.25, and logMi = 13.1, the best-ht curve in their Eig. 2b is well reproduced. 


3D POWER SPECTRA 

In the halo model computation of the cross-correlation power spectrum (PS), the 3D-PS is split in the one-halo (P^*) and two- 
halo (P^*) components with P = P^''+P^^. Eor a derivation of the P'^ and P^* discussed in the equations below, see Eomengo 
& Regis (2014). We remind that 5, denote 7 -ray astrophysical emitters (BL Lac, ESRQ, mAGN, and SEG), gj are associated to 
and galaxy populations (SDSS-DR 6 quasars, 2MASS galaxies, NVSS radio sources, SDSS-DR 8 Luminous Red Galaxies, and 
SDSS-DR 8 “main” galaxies), while 6 and stands for decaying and annihilating DM, respectively. In most of the equations, the 
dependence on z is not explicitely reported to simplify the notation. 

The 3D power spectrum of cross-correlation between 7 -rays from annihilating DM and galaxy catalogues is computed as: 


Pl'ldKz) 


P^max 


dM 


dn (Ng ) 
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Fig. 14.— Power spectrum (multiplied by k^) of cross-coiTelation between 7 -ray emitters and 2MASS galaxies at redshift z = 0.1. The left panel refers to 
annihilating DM. The right panel shows decaying DM and astrophysical sources. The matter power spectrum obtained within the halo model employed in this 
work is shown with a dashed black line and is compared with the halofit results (Takahashi et al. 2012) derived from high-resolution N-body simulations (black 
solid line). 




Fig. 15.— Power spectrum of cross-correlation between 7 -ray emitters and 2MASS galaxies at redshift z = 0.1. The power spectrum is multiplied by k and 
divided by the effective bias of 7 -ray emitters {b{z = 0.1)). The left panel compares predictions from annihilating (in the LOW scenario) and decaying DM with 
the astrophysical models of mAGN and SFG. The right panel focuses on mAGN and reports the power spectrum with three different assumptions for the M{C) 
relation (taken from Camera et al. (2014)). Thin lines show the 1-halo part of the power spectrum. 


The function u(k\M) is the Fourier transform of: 

u(x\M) = pl(x\M)/plj^ + bsub(M)ph(x\M)/M Jxpl(x\M)/ (C3) 

where p;, denotes the main halo profile and fejui, is the boost function associated to subhalos (introduced above). Note that 
u(k = 0\M) = Jd^xp^(x\M)/fP'. The product {Ngj)vg(k\m) is instead the Fourier transform of {NcsnjiM)) S^(x) + 

(Vsat,j(^)) Ph{x\M) / M . We have {Ng^)vg(k = Q\m) = (Ngj). The average number of galaxies gj at a given redshift is given by 
figjiz) = f dMdn/dM {Ngj). The details of the models of {Ngj), dn/dM and pi,(x\M) have been described in the previous Sections. 
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Fig. 16.— Left: Three dimensional power spectrum of cross-correlation between 7 -rays from annihilating DM (in the LOW scenario) and galaxies, evaluated 
at the redshift corresponding to the peak of the dNj/dz of each catalogue. Right: Same as left panel, but for SFG instead of DM. 


The impact of clustering assumptions on the 3D PS are illustrated in Fig. 14a, where we consider the 2MASS catalogue and 
show Pg. g 2 (k,z = 0.1). The boost from substructures makes the 7 -ray contributions from most massive halos to dominate the 
signal, and this is more pronounced in the HIGH case rather than in the LOW scenario. In the case without substructures, low mass 
halos becomes more important in the total budget of the 7 -ray emission. This explains the hierarchy at k ~ 1 /Mpc. For the same 
reasons, an opposite hierarchy occurs at very small scales (k > 100/Mpc). 

In the case of decaying DM, the PS of cross-correlation takes the form; 


Plj,siKz) = 
Pstsik,z) = 


dn {Ngj) 


/ dM 

dM n 

• rMmw. 


Vg(k\M)vs(k\M) 
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dn 
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(C4) 

(C5) 


Here ^, 5 (^ 144 ) is the Fourier transform of Ph{x\M )/In Fig. 14b, we show Pg. s{k,z = Q.l) (again for the 2MASS case), together 
with the matter power spectrum derived within our halo model approach. The latter is compared to a revised halofit PS derived 
from latest high-resolution N-body simulations (Takahashi et al. 2012). They agree within 20% at k < 10/Mpc and this supports 
our choices for the halo model ingredients. At larger k there is a departure, with less power in the halo model, but the picture at 
such small scales is in any case very uncertain, also from the simulations point of view. 

We assume astrophysical 7 -ray emitters to be point-like sources with the density field given by fs.{x-x') = Cg. 5^(x-x'). The 
3D PS of cross-correlation with galaxy catalogues can be written as: 
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(C 6 ) 

(C7) 


where bg, is the bias of 7 -ray astrophysical sources with respect to matter, for which we adopt bg.(C) = bf,{M{C)). Both Eqs. (C 6 ) 
and (C7) require the specification of the relation M(C) between the mass of the host halo M and the luminosity of the hosted 
object C. We will use the modeling of M(C) derived in Camera et al. (2014), where this aspect is discussed, and to which we 
refer the reader for the details. The blazar M(C) model of Camera et al. (2014) is adopted for both BL Lac and FSRQ. 

We caution that Eq. (C 6 ) for P^g gives only an approximate estimate of the 1-halo correlation. Indeed, modeling the satellite 
galaxies as a smooth component reduces their correlation with point-like 7 -ray sources. On the other hand, we assume that a 
halo hosting a given 7 -ray emitter also hosts the galaxies of all catalogues. This may not be true (e.g. some catalogue is mostly 
formed by galactic objects which do not host an AGN), thus artificially enhancing P^g.- Moreover, Eq. (C 6 ) is based on average 
relations, whilst a relative small number of outliers (i.e., bright 7 -ray sources in a halo with galaxies) can have a relevant impact. 
Eor all these reasons, and since P^g. is approximately independent on k, we can include in the ht an arbitrary constant term 

































22 


Cuoco et al. 



V 

1/3 





redshift z 


Fig. 17. — Effective bias for 7 -ray astrophysical emitters (left) and galaxies (right), as defined in Eqs. (C 8 ) and (C9), respectively. To illustrate the impact of 
the M(C) description, we additionally show the bias of mAGN when assuming the lower limit discussed in Camera et al. (2014) for such relation (green thin 
line). For comparison, we also report the bias of 7 -ray blazars considered in Xia et al. (2015) (red dotted line). In the galaxy cases (right panel), we show with 
circles the value of the different bias parameters adopted in Xia et al. (2015), where they were taken to be constant in redshift. The position of the dots refers to 
the redshift which corresponds to the peak of the galaxy distribution dNj/dz. 


allowing for both positive and negative corrections to Eq. (C 6 ). We call this additional quantity one-halo correction term, and we 
perform the analysis under the assumption that this term is not relevant (i.e. by setting it to vanish) and under the assumption that 
it is present, leaving it as a free parameter, one for each LSS tracer. 

Fig. 14 (right panel) shows Pg.,s,{k,z = 0.1), again taking the 2MASS catalogue as illustrative. The different classes of 7 -ray 
emitters show a similar spectrum, and have less (more) power than in the DM cases at intermediate (small) scales, as expected 
given their size. 

In Fig. 15, we show the difference between the DM and astrophysical PS at low redshift arising from the 1-halo term. To 
this aim we divide the PS by the bias in order to have the large scale PS (i.e., the two halo term) with a common normalization. 
At small scales the power associated to astrophysical sources is larger than for DM. The picture is opposite at intermediate 
scales, around Mpc, especially for SFG. The adopted model of M(C) makes the mAGN individual objects that contribute more 
to mAGN emission to be hosted in relatively large halos. This implies that the mAGN PS at Mpc scales is similar to the one of 
DM, explaining (part of) the origin of the degeneracy between mAGN and DM mentioned in the main text. We investigate the 
impact of different M{C) relations, taken from Camera et al. (2014), in the right panel of Fig. 15. 

The 3D PS of cross-correlation with the other catalogues are shown in Fig. 16. As illustrative examples, we selected the LOW 
scenario for annihilating DM, and SFG for astrophysical 7 -ray sources. The PS are computed at the redshift corresponding to the 
peak of the dNj/dz of each catalogue. 

In Fig. 17, we show the effective bias of astrophysical 7 -ray emitters and galaxies. They are defined with Pf/si ~ i^Si) 
at A: = 0, so they read; 


{bgjiz)) = 


if Si) 


pCmax,iiz) p 

/ dL<^i{ii,z)hm{cPj) 

min, 1 ( 2 ) 

dn (NgXM)) 

/ dM—h(M y 

dM ng. 


(C 8 ) 

(C9) 


Eq. (C 8 ) depends on the mass-luminosity relation M{C), while Eq. (C9) is governed by the modeling of {Ng.{M)) . The fair 
agreement shown by the computed bias with findings of autocorrelation studies quoted in the literature (e.g., Zehavi et al. (2005); 
Reid et al. (2014); Ho et al. (2012); White et al. (2011); Ross et al. (2010); Allevato et al. (2014)) is an important check of our 
modeling of M(C) and {Ng.(M)). 

The bias of 7 -ray blazars appears systematically lower than findings in Allevato et al. (2014). Translating the halo bias in terms 
of the mean mass hosting the blazars by means of Sheth cfe Tormen (1999), their results imply halos of M ~ 3 • lO'^M©, while 
according to our results shown in Fig. 17, FSRQs reside in halos of M ~ 1.5 • IO'^Mq and BL Lacs in M ~ 5 • lO'^M©. This is 
not surprising, if we consider that the work of Allevato et al. (2014) focuses on resolved objects, namely on a blazar subsample 
given by the brightest ones, which reside in more massive halos, while on the contrary, we investigate the unresolved component, 
which should be hosted by less massive halos. Moreover, the relatively low number of known 7 -ray objects prevents a firm 
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knowledge of their clustering, and sizable uncertainties on the bias are currently present. 

PLOTS OF CCF AT OTHER ENERGIES 
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Fig. 18.— Measured cross correlation function (CCF) (Xia et al. 2015) for E > 0.5 GeV, as a function of the angular separation 6 in the sky, compared to the 
best fit models of this analysis. The contribution to the CCF from the different astrophysical 7 -rays emitters (BL Lac, mAGN, SFG, FSRQ) are shown by dashed 
colored lines, while their sum (“Astro Total") and the DM contribution are indicated by solid green and red lines, respectively. The 1-halo correction term is 
shown as a solid blue line. The total contribution to the CCF is given by the black solid line. 
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Fig. 19.— Same as Fig. 18 but for £> 10 GeV. 
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